---
title: "Replication file: Divert when it does not hurt: The initiation of economic sanctions by US presidents from 1989 to 2015"
author: "Hana Attia"
date: "August 15, 2023"
output: html_document
---

```{r setup, include=FALSE}
rm(list = ls())
p_needed <- c("knitr", "MASS", "dplyr", "separationplot", "ggplot2", "reshape2", "ggthemes", "nnet",
              "mlogit", "Zelig", "stargazer", "haven", "foreign", "psych", "readr", "readxl", "haven", "Amelia", "sandwich", "lmtest", "sjPlot", "ggplot2", "sjmisc",               "patchwork", "haven")
packages <- rownames(installed.packages())
p_to_install <- p_needed[!(p_needed %in% packages)]
if (length(p_to_install) > 0) {
  install.packages(p_to_install)
}
sapply(p_needed, require, character.only = TRUE)
```

Load data
```{r}
#replace line below with your working directory
noICC_20210525 <- read_dta("/Users/hanaattia/Dropbox/Mac (2)/Desktop/PaperI_23042021/RoIE/R&R/RoIE_R&R_2/Submission_Replication/data_20221023.dta")
```

Figure 1: US executive sanctions by quarter (1989–2015)
```{r}
library(ggfortify)
plot1 <- noICC_20210525%>%
select(n_pres_imposition)%>%
ts(noICC_20210525$n_pres_imposition, frequency = 4, start = c(1989, 1), end = c(2015,4))
autoplot(plot1,
  xlab = "Year",
  ylab = "All sanctions decisions",
  ylim = c(0,8))

plot2 <- noICC_20210525%>%
select(n_pres_intense)%>%
ts(noICC_20210525$n_pres_intense, frequency = 4, start = c(1989, 1), end = c(2015,4))
autoplot(plot2,
  xlab = "Year",
  ylab = "Costly sanctions decisions",
  ylim = c(0,8))

plot3 <- noICC_20210525%>%
select(n_pres_not_intense)%>%
ts(noICC_20210525$n_pres_not_intense, frequency = 4, start = c(1989, 1), end = c(2015,4))
autoplot(plot3,
  xlab = "Year",
  ylab = "Non-costly sanctions decisions",
  ylim = c(0,8))
```

```{r}
noICC_20210525 <- noICC_20210525[order(noICC_20210525$quarter_year),]

noICC_20210525$lag_n_pres_imposition  <- lag(noICC_20210525$n_pres_imposition, n=1)
noICC_20210525$lag_n_pres_not_intense <- lag(noICC_20210525$n_pres_not_intense, n=1)
noICC_20210525$lag_n_pres_intense     <- lag(noICC_20210525$n_pres_intense, n=1)

noICC_20210525$lag_quarterly_unemp_unadjusted <- lag(noICC_20210525$quarterly_unemp_unadjusted, n=1)
noICC_20210525$lag_approval_ratings = lag(noICC_20210525$approval_ratings, n=1)
```

Table 2: Poisson models estimating the count of sanctions
```{r}
poi.sson1 <- glm(n_pres_imposition ~ lag_n_pres_imposition + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong +  lag_approval_ratings +  election_prox + EU_imposition, data = noICC_20210525, family = poisson)
summary(poi.sson1)

poi.sson2 <- glm(n_pres_not_intense ~ lag_n_pres_not_intense + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong +  lag_approval_ratings + election_prox + EU_imposition, data = noICC_20210525, family = poisson)
summary(poi.sson2)

poi.sson3 <- glm(n_pres_intense ~ lag_n_pres_intense + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong + lag_approval_ratings + election_prox + EU_imposition, data = noICC_20210525, family = poisson)
summary(poi.sson3)

poi.sson4 <- glm(n_pres_intense ~ lag_n_pres_intense + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong + lag_approval_ratings + election_prox + EU_intense, data = noICC_20210525, family = poisson)
summary(poi.sson4)
```

Figure 3: Marginal effect of interaction term on the imposition of non-costly sanctions
```{r}
poi.sson2 <- glm(n_pres_not_intense ~ lag_n_pres_not_intense + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong +  lag_approval_ratings + election_prox + EU_imposition, data = noICC_20210525, family = poisson)
summary(poi.sson2)


set_theme(
  base = theme_blank(),
  axis.title.size = .9,
  axis.textsize = .9,
  legend.size = .7,
  legend.title.size = .8
)

p <- plot_model(poi.sson2, type = "int", title = "", axis.title = c("Lagged unemployment (in %)", "Predicted count of non-costly sanctions"), legend.title = "President's party power", colors = "gs", grid = FALSE) + xlim(3.67, 10.5) + ylim(0,11)

p + scale_color_discrete(labels = c("Weak PPP", "Strong PPP"))  
```

Figure 4: Marginal effect of interaction term on the imposition of costly sanctions
```{r}
poi.sson4 <- glm(n_pres_intense ~ lag_n_pres_intense + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong + lag_approval_ratings + election_prox + EU_intense, data = noICC_20210525, family = poisson)
summary(poi.sson4)

set_theme(
  base = theme_blank(),
  axis.title.size = .9,
  axis.textsize = .9,
  legend.size = .7,
  legend.title.size = .8
)

p <- plot_model(poi.sson4, type = "int", title = "", axis.title = c("Logged nemployment (in %)", "Predicted count of costly sanctions"), legend.title = "President's party power", colors = "gs", grid = FALSE) + xlim(3.67, 10.5) + ylim(0,11)

p + scale_color_discrete(labels = c("Weak PPP", "Strong PPP"))  
```

#Appendix

Figure A1-1: Density plot of unemployment (in %, quarterly data)
```{r}
densityplot(noICC_20210525$quarterly_unemp_unadjusted, xlab = "Unemployment level (in %)", col = "palegreen4",
           ylab = "Kernel density")
```

Table A2-1: Negative Binomial models
```{r}
nb1 <- glm.nb(n_pres_imposition ~ lag_n_pres_imposition + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong +  lag_approval_ratings +  election_prox + EU_imposition, data = noICC_20210525)
summary(nb1)

nb2 <- glm.nb(n_pres_not_intense ~ lag_n_pres_not_intense + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong +  lag_approval_ratings + election_prox + EU_imposition, data = noICC_20210525)
summary(nb2)

nb3 <- glm.nb(n_pres_intense ~ lag_n_pres_intense + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong + lag_approval_ratings + election_prox + EU_imposition, data = noICC_20210525)
summary(nb3)

nb4 <- glm.nb(n_pres_intense ~ lag_n_pres_intense + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong + lag_approval_ratings + election_prox + EU_intense, data = noICC_20210525)
summary(nb4)
```

Figure A2-1.: Marginal effect of interaction term on the predicted count of all sanctions
```{r}
poi.sson1 <- glm(n_pres_imposition ~ lag_n_pres_imposition + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong +  lag_approval_ratings +  election_prox + EU_imposition, data = noICC_20210525, family = poisson)
summary(poi.sson1)

set_theme(
  base = theme_blank(),
  axis.title.size = .9,
  axis.textsize = .9,
  legend.size = .7,
  legend.title.size = .8
)

p <- plot_model(poi.sson1, type = "int", title = "", axis.title = c("Lagged unemployment (in %)", "Predicted count of all sanctions"), legend.title = "President's party power", colors = "gs", grid = FALSE) + xlim(3.67, 10.5) + ylim(0,11)

p + scale_color_discrete(labels = c("Weak PPP", "Strong PPP"))
```

Figure A2-2: Marginal effect of interaction term on the predicted count of non-costly sanctions
```{r}
poi.sson2 <- glm(n_pres_not_intense ~ lag_n_pres_not_intense + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong +  lag_approval_ratings + election_prox + EU_imposition, data = noICC_20210525, family = poisson)
summary(poi.sson2)

quantile(noICC_20210525$pres_party_power_cong, 0.25) #-11.87175
quantile(noICC_20210525$pres_party_power_cong, 0.75) #5.464836 

median(noICC_20210525$lag_n_pres_not_intense)
median(noICC_20210525$EU_imposition) #0
median(noICC_20210525$election_prox) #9
mean(noICC_20210525$lag_approval_ratings, na.rm = TRUE) #0.519101

par(mfrow=c(1,2))
newdat <- data.frame(lag_n_pres_not_intense = 1, lag_quarterly_unemp_unadjusted = 0:11, pres_party_power_cong = 5.464836, lag_approval_ratings=0.519101 , election_prox=9, EU_imposition = 0)
ginv <- poi.sson2$family$linkinv  ## inverse link function
prs <- predict(poi.sson2, newdata = newdat, type = "link", se.fit=TRUE)
newdat$pred <- ginv(prs[[1]])
newdat$lo <- ginv(prs[[1]] - 1.96 * prs[[2]])
newdat$up <- ginv(prs[[1]] + 1.96 * prs[[2]])
with(newdat, plot(lag_quarterly_unemp_unadjusted, pred, type = "l", ylim = c(min(lo), 6), xlab="Lagged unemployment (in %)", ylab= "Predicted count of non-costly sanctions", main = "Presidential party power high" ))
with(newdat, lines(lag_quarterly_unemp_unadjusted, lo, lty = 2, xlab="un"))
with(newdat, lines(lag_quarterly_unemp_unadjusted, up, lty = 2))

newdat <- data.frame(lag_n_pres_not_intense = 1, lag_quarterly_unemp_unadjusted = 0:11, pres_party_power_cong = -11.87175, lag_approval_ratings=0.519101 , election_prox=9, EU_imposition = 0)
ginv <- poi.sson2$family$linkinv  ## inverse link function
prs <- predict(poi.sson2, newdata = newdat,type = "link", se.fit=TRUE)
newdat$pred <- ginv(prs[[1]])
newdat$lo <- ginv(prs[[1]] - 1.96 * prs[[2]])
newdat$up <- ginv(prs[[1]] + 1.96 * prs[[2]])
with(newdat, plot(lag_quarterly_unemp_unadjusted, pred, type = "l", ylim = c(0, 6), xlab="Lagged unemployment (in %)", ylab= "Predicted count of non-costly sanctions", main = "Presidential party power low" ))
with(newdat, lines(lag_quarterly_unemp_unadjusted, lo, lty = 2))
with(newdat, lines(lag_quarterly_unemp_unadjusted, up, lty = 2))
```

A3.1 Autoregressive Poisson: see Stata for Table A3.1-1 and Figure A3.1-1

Table A3.2-1: Excluding sanctions against countries not prominent in the US public discourse (those below the 25th percentile)
```{r}
#replace line below with your working directory
noICC_20210525_25quantile <- read_dta("/Users/hanaattia/Dropbox/Mac (2)/Desktop/PaperI_23042021/RoIE/R&R/RoIE_R&R_2/Submission_Replication/data_20221023_NYT_25quant.dta")

noICC_20210525_25quantile <- noICC_20210525_25quantile[order(noICC_20210525_25quantile$quarter_year),]

noICC_20210525_25quantile$lag_n_pres_imposition  <- lag(noICC_20210525_25quantile$n_pres_imposition, n=1)
noICC_20210525_25quantile$lag_n_pres_not_intense <- lag(noICC_20210525_25quantile$n_pres_not_intense, n=1)
noICC_20210525_25quantile$lag_n_pres_intense     <- lag(noICC_20210525_25quantile$n_pres_intense, n=1)

noICC_20210525_25quantile$lag_quarterly_unemp_unadjusted <- lag(noICC_20210525_25quantile$quarterly_unemp_unadjusted, n=1)
noICC_20210525_25quantile$lag_approval_ratings = lag(noICC_20210525_25quantile$approval_ratings, n=1)

#Poisson
poi.sson1 <- glm(n_pres_imposition ~ lag_n_pres_imposition + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong +  lag_approval_ratings +  election_prox + EU_imposition, data = noICC_20210525_25quantile, family = poisson)
summary(poi.sson1)
poi.sson2 <- glm(n_pres_not_intense ~ lag_n_pres_not_intense + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong +  lag_approval_ratings + election_prox + EU_imposition, data = noICC_20210525_25quantile, family = poisson)
summary(poi.sson2)
poi.sson3 <- glm(n_pres_intense ~ lag_n_pres_intense + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong + lag_approval_ratings + election_prox + EU_intense, data = noICC_20210525_25quantile, family = poisson)
summary(poi.sson3)

#Negative binomial
#All sanctions
nb1 <- glm.nb(n_pres_imposition ~ lag_n_pres_imposition + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong +  lag_approval_ratings +  election_prox + EU_imposition, data = noICC_20210525_25quantile)
summary(nb1)
#Non-costly sanctions
nb2 <- glm.nb(n_pres_not_intense ~ lag_n_pres_not_intense + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong +  lag_approval_ratings + election_prox + EU_imposition, data = noICC_20210525_25quantile)
summary(nb2)
#Costly sanctions
nb3 <- glm.nb(n_pres_intense ~ lag_n_pres_intense + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong + lag_approval_ratings + election_prox + EU_intense, data = noICC_20210525_25quantile)
summary(nb3)
```


Figure A3.2-1: Marginal effect of interaction term using models that exclude countries not prominent in the US public discourse (those below 25th percentile)
```{r}
nb2 <- glm.nb(n_pres_not_intense ~ lag_n_pres_not_intense + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong +  lag_approval_ratings + election_prox + EU_imposition, data = noICC_20210525_25quantile)
summary(nb2)


set_theme(
  base = theme_blank(),
  axis.title.size = .9,
  axis.textsize = .9,
  legend.size = .7,
  legend.title.size = .8
)

p <- plot_model(nb2, type = "int", title = "", axis.title = c("Lagged unemployment (in %)", "Predicted count of non-costly sanctions"), legend.title = "President's party power", colors = "gs", grid = FALSE) + xlim(3.67, 10.5) + ylim(0,11)

p + scale_color_discrete(labels = c("Weak PPP", "Strong PPP"))  
```

Figure A3.2-2: Marginal effect of interaction term using models that exclude countries not prominent in the US public discourse (those below 25th percentile)
```{r}
nb2 <- glm.nb(n_pres_not_intense ~ lag_n_pres_not_intense + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong +  lag_approval_ratings + election_prox + EU_imposition, data = noICC_20210525_25quantile)
summary(nb2)


quantile(noICC_20210525_25quantile$pres_party_power_cong, 0.25) #-11.87175
quantile(noICC_20210525_25quantile$pres_party_power_cong, 0.75) #5.464836 

median(noICC_20210525_25quantile$lag_n_pres_not_intense)
median(noICC_20210525_25quantile$EU_imposition) #0
median(noICC_20210525_25quantile$election_prox) #9
mean(noICC_20210525_25quantile$lag_approval_ratings) #0.519101

par(mfrow=c(1,2))
newdat <- data.frame(lag_n_pres_not_intense = 1, lag_quarterly_unemp_unadjusted = 0:11, pres_party_power_cong = 5.464836, lag_approval_ratings=0.519101 , election_prox=9, EU_imposition = 0)
ginv <- nb2$family$linkinv  ## inverse link function
prs <- predict(nb2, newdata = newdat, type = "link", se.fit=TRUE)
newdat$pred <- ginv(prs[[1]])
newdat$lo <- ginv(prs[[1]] - 1.96 * prs[[2]])
newdat$up <- ginv(prs[[1]] + 1.96 * prs[[2]])
with(newdat, plot(lag_quarterly_unemp_unadjusted, pred, type = "l", ylim = c(min(lo), 6), xlab="Lagged unemployment (in %)", ylab= "Predicted count of non-costly sanctions", main = "Presidential party power high" ))
with(newdat, lines(lag_quarterly_unemp_unadjusted, lo, lty = 2, xlab="un"))
with(newdat, lines(lag_quarterly_unemp_unadjusted, up, lty = 2))

newdat <- data.frame(lag_n_pres_not_intense = 1, lag_quarterly_unemp_unadjusted = 0:11, pres_party_power_cong = -11.87175, lag_approval_ratings=0.519101 , election_prox=9, EU_imposition = 0)
ginv <- nb2$family$linkinv  ## inverse link function
prs <- predict(nb2, newdata = newdat,type = "link", se.fit=TRUE)
newdat$pred <- ginv(prs[[1]])
newdat$lo <- ginv(prs[[1]] - 1.96 * prs[[2]])
newdat$up <- ginv(prs[[1]] + 1.96 * prs[[2]])
with(newdat, plot(lag_quarterly_unemp_unadjusted, pred, type = "l", ylim = c(0, 6), xlab="Lagged unemployment (in %)", ylab= "Predicted count of non-costly sanctions", main = "Presidential party power low" ))
with(newdat, lines(lag_quarterly_unemp_unadjusted, lo, lty = 2))
with(newdat, lines(lag_quarterly_unemp_unadjusted, up, lty = 2))
```

Table A3.2-2: Excluding sanctions against countries not prominent in the US public discourse (those below the 50th percentile)
```{r}
#replace line below with your working directory
noICC_20210525_50quantile <- read_dta("/Users/hanaattia/Dropbox/Mac (2)/Desktop/PaperI_23042021/RoIE/R&R/RoIE_R&R_2/Submission_Replication/data_20221023_NYT_50quant.dta")

noICC_20210525_50quantile <- noICC_20210525_50quantile[order(noICC_20210525_50quantile$quarter_year),]

noICC_20210525_50quantile$lag_n_pres_imposition  <- lag(noICC_20210525_50quantile$n_pres_imposition, n=1)
noICC_20210525_50quantile$lag_n_pres_not_intense <- lag(noICC_20210525_50quantile$n_pres_not_intense, n=1)
noICC_20210525_50quantile$lag_n_pres_intense     <- lag(noICC_20210525_50quantile$n_pres_intense, n=1)

noICC_20210525_50quantile$lag_quarterly_unemp_unadjusted <- lag(noICC_20210525_50quantile$quarterly_unemp_unadjusted, n=1)
noICC_20210525_50quantile$lag_approval_ratings = lag(noICC_20210525_50quantile$approval_ratings, n=1)

poi.sson1 <- glm(n_pres_imposition ~ lag_n_pres_imposition + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong +  lag_approval_ratings +  election_prox + EU_imposition, data = noICC_20210525_50quantile, family = poisson)
summary(poi.sson1)

poi.sson2 <- glm(n_pres_not_intense ~ lag_n_pres_not_intense + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong +  lag_approval_ratings + election_prox + EU_imposition, data = noICC_20210525_50quantile, family = poisson)
summary(poi.sson2)

poi.sson3 <- glm(n_pres_intense ~ lag_n_pres_intense + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong + lag_approval_ratings + election_prox + EU_intense, data = noICC_20210525_50quantile, family = poisson)
summary(poi.sson3)
```

Figure A3.2-3: Marginal effect of interaction term using models that exclude countries not prominent in the US public discourse (those below 50th percentile)
```{r}
#install.packages('devtools', dependecies=TRUE)
#library(devtools)
#install_github("thomasp85/patchwork")

poi.sson4 <- glm(n_pres_not_intense ~ lag_n_pres_not_intense + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong +  lag_approval_ratings + election_prox + EU_imposition, data = noICC_20210525_50quantile, family = poisson)
summary(poi.sson4)


set_theme(
  base = theme_blank(),
  axis.title.size = .9,
  axis.textsize = .9,
  legend.size = .7,
  legend.title.size = .8
)

p <- plot_model(poi.sson4, type = "int", title = "", axis.title = c("Lagged unemployment (in %)", "Predicted count of non-costly sanctions"), legend.title = "President's party power", colors = "gs", grid = FALSE) + xlim(3.67, 10.5) + ylim(0,11)

p + scale_color_discrete(labels = c("Weak PPP", "Strong PPP"))  
```

Figure A3.2-4: Marginal effect of interaction term using models that exclude countries not prominent in the US public discourse (those below 50th percentile)
```{r}
poi.sson4 <- glm(n_pres_not_intense ~ lag_n_pres_not_intense + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong +  lag_approval_ratings + election_prox + EU_imposition, data = noICC_20210525_50quantile, family = poisson)
summary(poi.sson4)


quantile(noICC_20210525_50quantile$pres_party_power_cong, 0.25) #-11.87175
quantile(noICC_20210525_50quantile$pres_party_power_cong, 0.75) #5.464836 

sum((noICC_20210525_50quantile$n_pres_imposition))
median(noICC_20210525_50quantile$lag_n_pres_not_intense)
median(noICC_20210525_50quantile$EU_imposition) #0
median(noICC_20210525_50quantile$election_prox) #9
mean(noICC_20210525_50quantile$lag_approval_ratings) #0.519101

par(mfrow=c(1,2))
newdat <- data.frame(lag_n_pres_not_intense = 1, lag_quarterly_unemp_unadjusted = 0:11, pres_party_power_cong = 5.464836, lag_approval_ratings=0.519101 , election_prox=9, EU_imposition = 0)
ginv <- poi.sson4$family$linkinv  ## inverse link function
prs <- predict(poi.sson4, newdata = newdat, type = "link", se.fit=TRUE)
newdat$pred <- ginv(prs[[1]])
newdat$lo <- ginv(prs[[1]] - 1.96 * prs[[2]])
newdat$up <- ginv(prs[[1]] + 1.96 * prs[[2]])
with(newdat, plot(lag_quarterly_unemp_unadjusted, pred, type = "l", ylim = c(min(lo), 6), xlab="Lagged unemployment (in %)", ylab= "Predicted count of non-costly sanctions", main = "Presidential party power high" ))
with(newdat, lines(lag_quarterly_unemp_unadjusted, lo, lty = 2, xlab="un"))
with(newdat, lines(lag_quarterly_unemp_unadjusted, up, lty = 2))

newdat <- data.frame(lag_n_pres_not_intense = 1, lag_quarterly_unemp_unadjusted = 0:11, pres_party_power_cong = -11.87175, lag_approval_ratings=0.519101 , election_prox=9, EU_imposition = 0)
ginv <- poi.sson4$family$linkinv  ## inverse link function
prs <- predict(poi.sson4, newdata = newdat,type = "link", se.fit=TRUE)
newdat$pred <- ginv(prs[[1]])
newdat$lo <- ginv(prs[[1]] - 1.96 * prs[[2]])
newdat$up <- ginv(prs[[1]] + 1.96 * prs[[2]])
with(newdat, plot(lag_quarterly_unemp_unadjusted, pred, type = "l", ylim = c(0, 6), xlab="Lagged unemployment (in %)", ylab= "Predicted count of non-costly sanctions", main = "Presidential party power low" ))
with(newdat, lines(lag_quarterly_unemp_unadjusted, lo, lty = 2))
with(newdat, lines(lag_quarterly_unemp_unadjusted, up, lty = 2))
```

Table A3.3-1: Controlling for partisanship
```{r}
poi.sson1 <- glm(n_pres_imposition ~ lag_n_pres_imposition + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong +  lag_approval_ratings +  election_prox + EU_imposition + democrat, data = noICC_20210525, family = poisson)
summary(poi.sson1)

poi.sson2 <- glm(n_pres_not_intense ~ lag_n_pres_not_intense + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong +  lag_approval_ratings + election_prox + EU_imposition + democrat, data = noICC_20210525, family = poisson)
summary(poi.sson2)

poi.sson3 <- glm(n_pres_intense ~ lag_n_pres_intense + lag_quarterly_unemp_unadjusted + pres_party_power_cong + lag_quarterly_unemp_unadjusted*pres_party_power_cong + lag_approval_ratings + election_prox + EU_intense + democrat, data = noICC_20210525, family = poisson)
summary(poi.sson3)
```

A.3.4: Robustness check using clustered SE on the US administration level: See Stata for Table A3.4-1

A.3.5 Target dynamics: Loading data for dyadic analysis
```{r}
library(foreign)
rm(list = ls())

#replace line below with your working directory
data <- read_excel("/Users/hanaattia/Dropbox/Mac (2)/Desktop/PaperI_23042021/RoIE/R&R/RoIE_R&R_2/Submission_Replication/dyadic_analysis_20221106.xls")
#table(data$ccode)

data <- subset(data, !(data$ccode == 315 | data$ccode == 678 | data$ccode == 680))

data$t2_imposition <- data$t1_imposition^2
data$t3_imposition <- (data$t1_imposition^3)/1000

data$democrat <- 0 
data$democrat[data$year > 1992 & data$year < 2001] <- 1
data$democrat[data$year > 2008 & data$year < 2015] <- 1

data$bushsr <- 0
data$bushsr[data$year > 1989 & data$year < 1993] <- 1 

data$clinton <- 0
data$clinton[data$year > 1992 & data$year < 2001] <- 1

data$bushjr <- 0
data$bushjr[data$year > 2000 & data$year < 2009] <- 1

data$obama <- 0
data$obama[data$year > 2008 & data$year < 2015] <- 1
```

Figure A3.5-1: Density plot of unemployment (in %, annual data)
```{r}
densityplot(data$annual_unemp_unadjusted, xlab = "Unemployment level (in %)", col = "palegreen4",
           ylab = "Kernel density")
```

Table A3.5-1: Logit model
```{r}
library(miceadds)
model1 <- glm.cluster(data=data, formula=imposition ~  annual_unemp_unadjusted_lagged + pres_party_power_cong + annual_unemp_unadjusted_lagged*pres_party_power_cong + approval_lagged+ presidential_elections + logged_annual_trade_lagged + logged_eco_aid_lagged + logged_mil_aid_lagged + v2x_libdem_lagged + log_NYT_yearly_t_1 + t1_imposition + t2_imposition + t3_imposition, cluster="ccode", family="binomial")
summary(model1)

model2 <- glm.cluster(data=data, formula=imposition_not_intense ~  annual_unemp_unadjusted_lagged + pres_party_power_cong + annual_unemp_unadjusted_lagged*pres_party_power_cong + approval_lagged+ presidential_elections + logged_annual_trade_lagged + logged_eco_aid_lagged + logged_mil_aid_lagged + v2x_libdem_lagged + log_NYT_yearly_t_1 + t1_imposition + t2_imposition + t3_imposition, cluster="ccode", family="binomial")
summary(model2)

model3 <- glm.cluster(data=data, formula=imposition_intense ~  annual_unemp_unadjusted_lagged + pres_party_power_cong + annual_unemp_unadjusted_lagged*pres_party_power_cong + approval_lagged+ presidential_elections + logged_annual_trade_lagged + logged_eco_aid_lagged + logged_mil_aid_lagged + v2x_libdem_lagged + log_NYT_yearly_t_1 + t1_imposition + t2_imposition + t3_imposition, cluster="ccode", family="binomial")
summary(model3)
```

Figure A3.5-2: Non-costly sanctions, simulated marginal effect of interaction term
```{r}
data.sim.model4 <- subset(data, select = c(ccode, year, imposition_not_intense, annual_unemp_unadjusted_lagged, pres_party_power_cong, approval_lagged, presidential_elections,  logged_annual_trade_lagged, logged_eco_aid_lagged, logged_mil_aid_lagged, v2x_libdem_lagged, log_NYT_yearly_t_1, t1_imposition, t2_imposition, t3_imposition))
data.sim.model4 <-  data.sim.model4[complete.cases(data.sim.model4), ]

library(miceadds)
model2 <- glm.cluster(data=data.sim.model4, formula=imposition_not_intense ~  annual_unemp_unadjusted_lagged + pres_party_power_cong + annual_unemp_unadjusted_lagged*pres_party_power_cong + approval_lagged+ presidential_elections + logged_annual_trade_lagged + logged_eco_aid_lagged + logged_mil_aid_lagged + v2x_libdem_lagged + log_NYT_yearly_t_1 + t1_imposition + t2_imposition + t3_imposition, cluster="ccode", family="binomial")
summary(model2)

gamma.hat2 <- coef(model2)
V.hat2 <- vcov(model2)

set.seed(345)
library(MASS)
S2 <- mvrnorm(1000,gamma.hat2,V.hat2)

unemp.min <- min(data.sim.model4$annual_unemp_unadjusted_lagged)
unemp.q1 <- quantile(data.sim.model4$annual_unemp_unadjusted_lagged, 0.25)
unemp.mean <- mean(data.sim.model4$annual_unemp_unadjusted_lagged)
unemp.q3 <- quantile(data.sim.model4$annual_unemp_unadjusted_lagged, 0.75)
unemp.max <- max(data.sim.model4$annual_unemp_unadjusted_lagged)

ppp.min <- min(data.sim.model4$pres_party_power_cong)
ppp.q1 <- quantile(data.sim.model4$pres_party_power_cong, 0.25)
ppp.max <- max(data.sim.model4$pres_party_power_cong)
ppp.q3 <- quantile(data.sim.model4$pres_party_power_cong, 0.75)

sim.unemp <-seq(min(data.sim.model4$annual_unemp_unadjusted_lagged), max(data.sim.model4$annual_unemp_unadjusted_lagged),length.out = 100)

scenario.pppmin<- cbind(1, sim.unemp, ppp.min, mean(data.sim.model4$approval_lagged), mean(data.sim.model4$presidential_elections), mean(data.sim.model4$logged_annual_trade_lagged), mean(data.sim.model4$logged_eco_aid_lagged), mean(data.sim.model4$logged_mil_aid_lagged), mean(data.sim.model4$v2x_libdem_lagged), mean(data.sim.model4$log_NYT_yearly_t_1), median(data.sim.model4$t1_imposition), median(data.sim.model4$t2_imposition), median(data.sim.model4$t3_imposition), sim.unemp*ppp.min)
scenario.pppmax <- cbind(1, sim.unemp, ppp.max,mean(data.sim.model4$approval_lagged), mean(data.sim.model4$presidential_elections), mean(data.sim.model4$logged_annual_trade_lagged), mean(data.sim.model4$logged_eco_aid_lagged), mean(data.sim.model4$logged_mil_aid_lagged), mean(data.sim.model4$v2x_libdem_lagged), mean(data.sim.model4$log_NYT_yearly_t_1), median(data.sim.model4$t1_imposition), median(data.sim.model4$t2_imposition), median(data.sim.model4$t3_imposition), sim.unemp*ppp.max)

Xbeta_pppmin <- S2 %*%(t(scenario.pppmin))
Xbeta_pppmax <- S2 %*%(t(scenario.pppmax))
p.sim.min <- (exp(Xbeta_pppmin))/ (1 + exp(Xbeta_pppmin))
p.sim.max <- (exp(Xbeta_pppmax))/ (1 + exp(Xbeta_pppmax))

fd1 <- p.sim.max - p.sim.min
p.mean.fd1 <- apply(fd1, 2, mean)
p.qu.fd1 <- t(apply(fd1, 2, quantile, prob = c(0.025, 0.975)))
p.mean.min <- apply(p.sim.min, 2, mean)
p.qu.min <- t(apply(p.sim.min, 2, quantile, prob = c(0.025, 0.975)))
p.mean.max <- apply(p.sim.max, 2, mean)
p.qu.max <- t(apply(p.sim.max, 2, quantile, prob = c(0.025, 0.975)))

results <- as.data.frame(p.mean.fd1)
results$free <- sim.unemp
results$lower <- p.qu.fd1 [,1]
results$upper <- p.qu.fd1 [,2]

plot((sim.unemp), p.mean.min, type="n",
ylim = c(0, 0.3),
cex=0.9,
xlab = "Lagged unemployment (in %)",
ylab = "Predicted probability of non-costly sanctions")
polygon(c(rev((sim.unemp)), (sim.unemp)), c(rev(p.qu.min[,2]), p.qu.min[,1]),
col = "lightblue",
border = NA)
polygon(c(rev((sim.unemp)), (sim.unemp)), c(rev(p.qu.max[,2]), p.qu.max[,1]),
col = "lightcoral",
border = NA)

lines((sim.unemp), p.mean.min, lwd = 1)
#lines((sim.unemp), p.qu.min[, 1], lty = "dashed", col = "gray20")
#lines((sim.unemp), p.qu.min[, 2], lty = "dashed", col = "gray20")
lines((sim.unemp), p.mean.max, lwd = 1)
#lines((sim.unemp), p.qu.max[, 1], lty = "dashed", col = "gray20")
#lines((sim.unemp), p.qu.max[, 2], lty = "dashed", col = "gray20")
legend("topright",
title="Presidential party power",
c("Low", "High"),
cex=0.8,
bty = 'n',
lty=c(1,1),
col=c("lightblue","lightcoral"),
xpd = TRUE)
```

Figure A3.5-3: Costly sanctions, simulated marginal effect of interaction term
```{r}
data.sim.model4 <- subset(data, select = c(ccode, year, imposition_intense, annual_unemp_unadjusted_lagged, pres_party_power_cong, approval_lagged, presidential_elections,  logged_annual_trade_lagged, logged_eco_aid_lagged, logged_mil_aid_lagged, v2x_libdem_lagged, log_NYT_yearly_t_1, t1_imposition, t2_imposition, t3_imposition))
data.sim.model4 <-  data.sim.model4[complete.cases(data.sim.model4), ]

library(miceadds)
model3 <- glm.cluster(data=data.sim.model4, formula=imposition_intense ~  annual_unemp_unadjusted_lagged + pres_party_power_cong + annual_unemp_unadjusted_lagged*pres_party_power_cong + approval_lagged+ presidential_elections + logged_annual_trade_lagged + logged_eco_aid_lagged + logged_mil_aid_lagged + v2x_libdem_lagged + log_NYT_yearly_t_1 + t1_imposition + t2_imposition + t3_imposition, cluster="ccode", family="binomial")
summary(model3)

gamma.hat2 <- coef(model3)
V.hat2 <- vcov(model3)

set.seed(345)
library(MASS)
S2 <- mvrnorm(1000,gamma.hat2,V.hat2)

unemp.min <- min(data.sim.model4$annual_unemp_unadjusted_lagged)
unemp.q1 <- quantile(data.sim.model4$annual_unemp_unadjusted_lagged, 0.25)
unemp.mean <- mean(data.sim.model4$annual_unemp_unadjusted_lagged)
unemp.q3 <- quantile(data.sim.model4$annual_unemp_unadjusted_lagged, 0.75)
unemp.max <- max(data.sim.model4$annual_unemp_unadjusted_lagged)

ppp.min <- min(data.sim.model4$pres_party_power_cong)
ppp.q1 <- quantile(data.sim.model4$pres_party_power_cong, 0.25)
ppp.max <- max(data.sim.model4$pres_party_power_cong)
ppp.q3 <- quantile(data.sim.model4$pres_party_power_cong, 0.75)

sim.unemp <-seq(min(data.sim.model4$annual_unemp_unadjusted_lagged), max(data.sim.model4$annual_unemp_unadjusted_lagged),length.out = 100)

scenario.pppmin<- cbind(1, sim.unemp, ppp.min, mean(data.sim.model4$approval_lagged), mean(data.sim.model4$presidential_elections), mean(data.sim.model4$logged_annual_trade_lagged), mean(data.sim.model4$logged_eco_aid_lagged), mean(data.sim.model4$logged_mil_aid_lagged), median(data.sim.model4$v2x_libdem_lagged), mean(data.sim.model4$log_NYT_yearly_t_1), median(data.sim.model4$t1_imposition), median(data.sim.model4$t2_imposition), median(data.sim.model4$t3_imposition), sim.unemp*ppp.min)
scenario.pppmax <- cbind(1, sim.unemp, ppp.max,mean(data.sim.model4$approval_lagged), mean(data.sim.model4$presidential_elections), mean(data.sim.model4$logged_annual_trade_lagged), mean(data.sim.model4$logged_eco_aid_lagged), mean(data.sim.model4$logged_mil_aid_lagged), median(data.sim.model4$v2x_libdem_lagged), mean(data.sim.model4$log_NYT_yearly_t_1), median(data.sim.model4$t1_imposition), median(data.sim.model4$t2_imposition), median(data.sim.model4$t3_imposition), sim.unemp*ppp.max)

Xbeta_pppmin <- S2 %*%(t(scenario.pppmin))
Xbeta_pppmax <- S2 %*%(t(scenario.pppmax))
p.sim.min <- (exp(Xbeta_pppmin))/ (1 + exp(Xbeta_pppmin))
p.sim.max <- (exp(Xbeta_pppmax))/ (1 + exp(Xbeta_pppmax))

fd1 <- p.sim.max - p.sim.min
p.mean.fd1 <- apply(fd1, 2, mean)
p.qu.fd1 <- t(apply(fd1, 2, quantile, prob = c(0.025, 0.975)))
p.mean.min <- apply(p.sim.min, 2, mean)
p.qu.min <- t(apply(p.sim.min, 2, quantile, prob = c(0.025, 0.975)))
p.mean.max <- apply(p.sim.max, 2, mean)
p.qu.max <- t(apply(p.sim.max, 2, quantile, prob = c(0.025, 0.975)))

results <- as.data.frame(p.mean.fd1)
results$free <- sim.unemp
results$lower <- p.qu.fd1 [,1]
results$upper <- p.qu.fd1 [,2]

plot((sim.unemp), p.mean.min, type="n",
ylim = c(0, 0.3),
cex=0.9,
xlab = "Lagged unemployment (in %)",
ylab = "Predicted probability of costly sanctions")
polygon(c(rev((sim.unemp)), (sim.unemp)), c(rev(p.qu.min[,2]), p.qu.min[,1]),
col = "lightblue",
border = NA)
polygon(c(rev((sim.unemp)), (sim.unemp)), c(rev(p.qu.max[,2]), p.qu.max[,1]),
col = "lightcoral",
border = NA)

lines((sim.unemp), p.mean.min, lwd = 1)
#lines((sim.unemp), p.qu.min[, 1], lty = "dashed", col = "gray20")
#lines((sim.unemp), p.qu.min[, 2], lty = "dashed", col = "gray20")
lines((sim.unemp), p.mean.max, lwd = 1)
#lines((sim.unemp), p.qu.max[, 1], lty = "dashed", col = "gray20")
#lines((sim.unemp), p.qu.max[, 2], lty = "dashed", col = "gray20")
legend("topright",
title="Presidential party power",
c("Low", "High"),
cex=0.8,
bty = 'n',
lty=c(1,1),
col=c("lightblue","lightcoral"),
xpd = TRUE)
```

Table A3.5-2: Logit model clustering SE on the presidential administration level
```{r}
library(miceadds)
model1 <- glm.cluster(data=data, formula=imposition ~  annual_unemp_unadjusted_lagged + pres_party_power_cong + annual_unemp_unadjusted_lagged*pres_party_power_cong + approval_lagged+ presidential_elections + logged_annual_trade_lagged + logged_eco_aid_lagged + logged_mil_aid_lagged + v2x_libdem_lagged + log_NYT_yearly_t_1 + t1_imposition + t2_imposition + t3_imposition, cluster="pres", family="binomial")
summary(model1)

model2 <- glm.cluster(data=data, formula=imposition_not_intense ~  annual_unemp_unadjusted_lagged + pres_party_power_cong + annual_unemp_unadjusted_lagged*pres_party_power_cong + approval_lagged+ presidential_elections + logged_annual_trade_lagged + logged_eco_aid_lagged + logged_mil_aid_lagged + v2x_libdem_lagged + log_NYT_yearly_t_1 +  t1_imposition + t2_imposition + t3_imposition, cluster="pres", family="binomial")
summary(model2)

model3 <- glm.cluster(data=data, formula=imposition_intense ~  annual_unemp_unadjusted_lagged + pres_party_power_cong + annual_unemp_unadjusted_lagged*pres_party_power_cong + approval_lagged+ presidential_elections + logged_annual_trade_lagged + logged_eco_aid_lagged + logged_mil_aid_lagged + v2x_libdem_lagged + log_NYT_yearly_t_1 + t1_imposition + t2_imposition + t3_imposition, cluster="pres", family="binomial")
summary(model3)
```

Table A3.5-3: Logit model controlling for EU sanctions
```{r}
library(miceadds)

model1 <- glm.cluster(data=data, formula=imposition ~  annual_unemp_unadjusted_lagged + pres_party_power_cong + annual_unemp_unadjusted_lagged*pres_party_power_cong + approval_lagged+ presidential_elections + logged_annual_trade_lagged + logged_eco_aid_lagged + logged_mil_aid_lagged + v2x_libdem_lagged + log_NYT_yearly_t_1 + EU_imposition +  t1_imposition + t2_imposition + t3_imposition, cluster="ccode", family="binomial")
summary(model1)

model2 <- glm.cluster(data=data, formula=imposition_not_intense ~  annual_unemp_unadjusted_lagged + pres_party_power_cong + annual_unemp_unadjusted_lagged*pres_party_power_cong + approval_lagged+ presidential_elections + logged_annual_trade_lagged + logged_eco_aid_lagged + logged_mil_aid_lagged + v2x_libdem_lagged + log_NYT_yearly_t_1 + EU_imposition + t1_imposition + t2_imposition + t3_imposition, cluster="ccode", family="binomial")
summary(model2)

model3 <- glm.cluster(data=data, formula=imposition_intense ~  annual_unemp_unadjusted_lagged + pres_party_power_cong + annual_unemp_unadjusted_lagged*pres_party_power_cong + approval_lagged+ presidential_elections + logged_annual_trade_lagged + logged_eco_aid_lagged + logged_mil_aid_lagged + v2x_libdem_lagged + log_NYT_yearly_t_1 + EU_imposition + t1_imposition + t2_imposition + t3_imposition, cluster="ccode", family="binomial")
summary(model3)
```
